%%author:sousanone
%%email:nj2000515@foxmail.com
%%date:2023/3/23



clc,clear all;

addpath 'D:\matlab2022\bin\cmi'

% addpath('D:\mat\bin\ptychography-master\utils')
%% FZA generation
lamda = 1030e-6;                                        % 入射波长，单位mm米、, 632.8nm nm-um-mm-m
k0=2*pi/lamda;                                        
Nx =600; 
range1=2000;
focallength=30;                                         %焦距 
step1=range1/(Nx-1);                                     % 步距
yp1=-range1/2:step1:range1/2;
xp1=-range1/2:step1:range1/2;
[x1,y1]=meshgrid(xp1,yp1);                               % 划分网格
position1  = x1 + 1i * y1;                               % 某点的位置
lensletr = abs(position1);                               % 某点的半径
lensang  = angle(position1);
Nc=floor(Nx/2)+1;
TFZP= exp(-1i*pi*lensletr.^2/(lamda*focallength));
a=sign(real(TFZP));
lens1 =sign(1+sign(real(TFZP)));   



%%
lamda = 1030e-6;
pix = 5.5e-3;
z1 = 20;
z2 = 30;
mask1=zeros(Nx,Nx);
mask1(Nx/2-200:Nx/2+200,Nx/2-200:Nx/2+200)=1;

% mask1=circle_mask(Nx*2,400,Nx/2,Nx/2);


% lens2=mask1.*lens1;
% % % 这啥也不是
% cmi_prior = lens2;%.*exp(1i*lens2/max(max(lens2)));

% % 相位型
cmi_prior = exp(1i*pi/8*lens1/max(max(lens1)));%.*mask1;

% % 振幅形
% cmi_prior = lens2/max(max(lens2)).*exp(1i*mask1);


obj=imread('cameraman.tif');
obj=imresize(obj,[Nx*2,Nx*2]);
obj1=double(obj);
phase=imread('mla.jpg');
phase=double(rgb2gray(phase));
phase=imresize(phase,[Nx*2,Nx*2]);
phase=phase;

% obj2=double(obj1).*exp(1i*phase/max(max(phase)));

obj2=phase.*exp(1i.*obj1/max(max(obj1)));




%% rawdata generation
s=double(imread('Lena.tif'));
s=imresize(s,[600,600]);
s=exp(1i.*s/max(max(s)));
cmi_prior=s;
probe= f_Gaussian_beams(1e-3,0,lamda,2e-3,2e-3,600,600);
sstep=150;
col=5;
row=5;
raw=zeros(600,600,row*col);
numd=1;
for i=1:row
    for j=1:col
        obj22=obj2(1+(i-1)*sstep:(i-1)*sstep+600,1+(j-1)*sstep:(j-1)*sstep+600).*probe;
        temp1=Propagate(obj22,'fresnel',pix,lamda,z1);
        temp2=temp1.*double(cmi_prior);
        temp3=Propagate(temp2,'fresnel',pix,lamda,z2);
        raw(:,:,numd)=abs(temp3);
%         imshow(raw(:,:,numd),[])
        numd=numd+1;
        
    end
end





%% CMI reconstruction
iter = 2;
obj_g = ones(Nx*2,2*Nx);
% probe22=zeros(Nx,Nx);
probe22 =abs(probe);
cmi_prior=ones(600,600);

for i = 1:iter
    for i1=1:row
        for j1=1:col

            outputimg1=obj_g(1+(i1-1)*sstep:(i1-1)*sstep+600,1+(j1-1)*sstep:(j1-1)*sstep+600);
%             imshow(outputimg1,[])
            ew1 = outputimg1.*probe22;
%             imshow(ew1,[])

            iw2 = Propagate(ew1,'fresnel',pix,lamda,z1);
%             imshow(iw2,[])

            ew2 = iw2.*double(cmi_prior);
%             imshow(ew2,[])

            dp1 = Propagate(ew2,'fresnel',pix,lamda,z2);
%             imshow(dp1,[])

            dp2 = raw(:,:,(i1-1)*row+j1).*exp(1i.*angle(dp1));
%             imshow(dp2,[])
    
            new2 = Propagate(dp2,'fresnel',pix,lamda,-z2);
%             imshow(new2,[])

            niw2 = iw2 + conj(double(cmi_prior))./(max(max(abs(double(cmi_prior)).^2))).*(new2-ew2);
%             imshow(niw2,[])

            cmi_prior=exp(1i.*angle(cmi_prior+(conj(iw2))./(max(max(abs(double(iw2)).^2))).*(new2-ew2)));
%             cmi_prior=cmi_prior.*mask1;
%             imshow(abs(cmi_prior),[])
 
            new1 = Propagate(niw2,'fresnel',pix,lamda,-z1);
%             imshow(new1,[])

            outputimg2= outputimg1+conj(probe22)./max(max(abs(probe22).^2)).*(new1-ew1);
%             imshow(outputimg2,[])

            probe22=probe22+conj(outputimg1)./max(max(abs(outputimg1).^2)).*(new1-ew1);
%             imshow(probe22,[])

            obj_g(1+(i1-1)*sstep:(i1-1)*sstep+600,1+(j1-1)*sstep:(j1-1)*sstep+600)=outputimg2;
%             imshow(abs(obj_g),[])
        end
    end
end

            


















subplot(321)
imshow(abs(obj_g),[])
title('物体振幅')
subplot(322)
imshow(angle(obj_g),[])
title('物体相位')
subplot(323)
imshow(abs(cmi_prior),[])
title('掩膜振幅')
subplot(324)
imshow(angle(cmi_prior),[])
title('掩膜相位')
subplot(325)
imshow(abs(probe22),[])
title('探针振幅')
subplot(326)
imshow(angle(probe22),[])
title('探针相位')
% % subplot(321)
% % imshow(abs(obj_g))
% % title('物体振幅')
% % subplot(322)
% % imshow(angle(obj_g))
% % title('物体相位')
% % subplot(323)
% % imshow(abs(cmi_prior))
% % title('掩膜振幅')
% % subplot(324)
% % imshow(angle(cmi_prior))
% % title('掩膜相位')
% % subplot(325)
% % imshow(abs(probe22))
% % title('探针振幅')
% % subplot(326)
% % imshow(angle(probe22))
% % title('探针相位')






